#include "VGQAsympExp.h"
#include "PolynomialTerms.h"
#include "Polynomial.h"
#include <math.h>
#include <boost/math/special_functions/gamma.hpp>
#include <iostream>

VGQAsympExp::VGQAsympExp(double mu, double lambda, double alpha, double beta, size_t nTerms) : _mu(mu), _lambda(lambda), _alpha(alpha), _beta(beta), _nTerms(nTerms)
{
	using namespace VGQuantile;
	using namespace boost;

	_const1 = pow(2*alpha,lambda) * pow((alpha*alpha)-(beta*beta),-lambda) * math::tgamma(lambda);


	// Construct Polynomials (u->1)
	//P_0
	Polynomial rp0;
	rp0.setCoefficient(0, c0_0(lambda, alpha, beta));
	rp0.setCoefficient(1, c0_1(lambda, alpha, beta));
	_rightPolynomials.push_back(rp0);

	//P_1
	Polynomial rp1;
	rp1.setCoefficient(0, c1_0(lambda, alpha, beta));
	rp1.setCoefficient(1, c1_1(lambda, alpha, beta));
	_rightPolynomials.push_back(rp1);

	//P_2
	Polynomial rp2;
	rp2.setCoefficient(0, c2_0(lambda, alpha, beta));
	rp2.setCoefficient(1, c2_1(lambda, alpha, beta));
	rp2.setCoefficient(2, c2_2(lambda, alpha, beta));
	_rightPolynomials.push_back(rp2);

	//P_3
	Polynomial rp3;
	rp3.setCoefficient(0, c3_0(lambda, alpha, beta));
	rp3.setCoefficient(1, c3_1(lambda, alpha, beta));
	rp3.setCoefficient(2, c3_2(lambda, alpha, beta));
	rp3.setCoefficient(3, c3_3(lambda, alpha, beta));
	_rightPolynomials.push_back(rp3);

	//P_4
	Polynomial rp4;
	rp4.setCoefficient(0, c4_0(lambda, alpha, beta));
	rp4.setCoefficient(1, c4_1(lambda, alpha, beta));
	rp4.setCoefficient(2, c4_2(lambda, alpha, beta));
	rp4.setCoefficient(3, c4_3(lambda, alpha, beta));
	rp4.setCoefficient(4, c4_4(lambda, alpha, beta));
	_rightPolynomials.push_back(rp4);

	//P_5
	Polynomial rp5;
	rp5.setCoefficient(0, c5_0(lambda, alpha, beta));
	rp5.setCoefficient(1, c5_1(lambda, alpha, beta));
	rp5.setCoefficient(2, c5_2(lambda, alpha, beta));
	rp5.setCoefficient(3, c5_3(lambda, alpha, beta));
	rp5.setCoefficient(4, c5_4(lambda, alpha, beta));
	rp5.setCoefficient(5, c5_5(lambda, alpha, beta));
	_rightPolynomials.push_back(rp5);
	
	// Construct Polynomials (u -> 0)
	//P_0
	Polynomial lp0;
	lp0.setCoefficient(0, c0_0(lambda, alpha, -beta));
	lp0.setCoefficient(1, c0_1(lambda, alpha, -beta));
	_leftPolynomials.push_back(lp0);

	//P_1
	Polynomial lp1;
	lp1.setCoefficient(0, c1_0(lambda, alpha, -beta));
	lp1.setCoefficient(1, c1_1(lambda, alpha, -beta));
	_leftPolynomials.push_back(lp1);

	//P_2
	Polynomial lp2;
	lp2.setCoefficient(0, c2_0(lambda, alpha, -beta));
	lp2.setCoefficient(1, c2_1(lambda, alpha, -beta));
	lp2.setCoefficient(2, c2_2(lambda, alpha, -beta));
	_leftPolynomials.push_back(lp2);

	//P_3
	Polynomial lp3;
	lp3.setCoefficient(0, c3_0(lambda, alpha, -beta));
	lp3.setCoefficient(1, c3_1(lambda, alpha, -beta));
	lp3.setCoefficient(2, c3_2(lambda, alpha, -beta));
	lp3.setCoefficient(3, c3_3(lambda, alpha, -beta));
	_leftPolynomials.push_back(lp3);

	//P_4
	Polynomial lp4;
	lp4.setCoefficient(0, c4_0(lambda, alpha, -beta));
	lp4.setCoefficient(1, c4_1(lambda, alpha, -beta));
	lp4.setCoefficient(2, c4_2(lambda, alpha, -beta));
	lp4.setCoefficient(3, c4_3(lambda, alpha, -beta));
	lp4.setCoefficient(4, c4_4(lambda, alpha, -beta));
	_leftPolynomials.push_back(lp4);

	//P_5
	Polynomial lp5;
	lp5.setCoefficient(0, c5_0(lambda, alpha, -beta));
	lp5.setCoefficient(1, c5_1(lambda, alpha, -beta));
	lp5.setCoefficient(2, c5_2(lambda, alpha, -beta));
	lp5.setCoefficient(3, c5_3(lambda, alpha, -beta));
	lp5.setCoefficient(4, c5_4(lambda, alpha, -beta));
	lp5.setCoefficient(5, c5_5(lambda, alpha, -beta));
	_leftPolynomials.push_back(lp5);

}


VGQAsympExp::~VGQAsympExp(void)
{
}

double VGQAsympExp::LeftTail( double u ) const
{
	return _mu - Evaluate(_leftPolynomials, -_beta, 1-u);	
}

double VGQAsympExp::RightTail( double u ) const
{
	return _mu + Evaluate(_rightPolynomials, _beta, u);	
}

double VGQAsympExp::yVar(double beta, double u) const
{
	double v = (1-u)*_const1;
	return -log(v)/(_alpha - beta);
}

double VGQAsympExp::Evaluate( const std::vector<Polynomial>& polynomials, double beta, double u ) const
{
	//TODO: Optimally truncate asymptotic expansion
	
	double y = yVar(beta, u);
	double z = log(y);
	
	double sum = y;
	for(int i = 0; i <= _nTerms; i++)
	{
		sum += ( polynomials[i].evaluate(z) / pow(y, i) );
	}

	return sum;
}
